COVID19 UK

Prepare data

Read and format data

Prevalence

df_uk_prev <- read_csv('UK_timeseries_prep_2005.csv')

Parsed with column specification:
cols(
  ut_area = col_character(),
  date = col_character(),
  cumcase = col_double(),
  poptotal = col_double(),
  rate = col_double()
)

df_uk_prev <- df_uk_prev %>% 
  select(ut_area, date, rate) %>% 
  rename(rate_day = rate) %>%
  mutate(date = as.Date(date, "%d%b%Y"))

df_uk_prev

Personality

df_uk_pers <- read_csv('timeseries_uk_utla_march9_april_09.csv')

Parsed with column specification:
cols(
  ut_area = col_character(),
  time = col_double(),
  areaname = col_character(),
  open = col_double(),
  extra = col_double(),
  agree = col_double(),
  neuro = col_double(),
  sci = col_double(),
  frequ = col_double(),
  ut_name = col_character(),
  poptotal = col_double(),
  rate_day = col_double()
)

df_uk_pers <- df_uk_pers %>% 
  select(ut_area, open, agree, neuro, sci, extra) %>% 
  dplyr::rename(pers_o = open, 
         pers_c = sci,
         pers_e = extra,
         pers_a = agree,
         pers_n = neuro) %>%
  distinct()

df_uk_pers

NA

Social distancing

df_uk_socdist <- read_csv('UK_socdist_fb_nuts3.csv')

Parsed with column specification:
cols(
  nuts3 = col_character(),
  date = col_date(format = ""),
  all_day_bing_tiles_visited_relat = col_double(),
  all_day_ratio_single_tile_users = col_double(),
  open = col_double(),
  extra = col_double(),
  agree = col_double(),
  neuro = col_double(),
  sci = col_double(),
  frequ = col_double(),
  nuts3_name = col_character(),
  runday = col_double()
)

df_uk_socdist$date %>% summary()

        Min.      1st Qu.       Median         Mean      3rd Qu.         Max. 
"2020-03-01" "2020-03-08" "2020-03-16" "2020-03-16" "2020-03-24" "2020-03-31" 

df_uk_socdist <- df_uk_socdist %>% select(-runday, -frequ) %>%
  dplyr::rename(pers_o = open, 
                pers_c = sci,
                pers_e = extra,
                pers_a = agree,
                pers_n = neuro) %>% 
  select(-nuts3_name) %>% 
  dplyr::rename(socdist_tiles = all_day_bing_tiles_visited_relat,
                socdist_single_tile = all_day_ratio_single_tile_users) %>%
  drop_na()

df_uk_socdist

Controls

df_uk_ctrl_nuts <- read_csv("controls_UK_nuts3.csv")

Parsed with column specification:
cols(
  nuts3 = col_character(),
  nuts3_name = col_character(),
  airport_dist = col_double(),
  males = col_double(),
  popdens = col_double(),
  manufacturing = col_double(),
  tourism = col_double(),
  health = col_double(),
  academic = col_double(),
  medinc = col_double(),
  medage = col_double(),
  conservative = col_double()
)

df_uk_ctrl_nuts <- df_uk_ctrl_nuts %>% select(-nuts3_name)
df_uk_ctrl_nuts

df_uk_ctrl_ut <- read_csv("controls_UK_ut.csv")

Parsed with column specification:
cols(
  ut_area = col_character(),
  ut_name = col_character(),
  airport_dist = col_double(),
  males = col_double(),
  popdens = col_double(),
  manufacturing = col_double(),
  tourism = col_double(),
  health = col_double(),
  academic = col_double(),
  medinc = col_double(),
  medage = col_double(),
  conservative = col_double()
)

df_uk_ctrl_ut <- df_uk_ctrl_ut %>% select(-ut_name)
df_uk_ctrl_ut

NA
NA

Merge prevalence data

df_uk <- df_uk_prev %>% 
  plyr::join(df_uk_pers, by='ut_area') %>% 
  plyr::join(df_uk_ctrl_ut, by='ut_area')

# create sequence of dates
date_sequence <- seq.Date(min(df_uk$date),
                          max(df_uk$date), 1)
                     
# create data frame with time sequence
df_dates = tibble(date_sequence, 1:length(date_sequence)) 
names(df_dates) <- c('date', 'time')

# merge day index with gps data
df_uk = df_uk %>% 
  merge(df_dates, by='date') %>% 
  arrange(ut_area) %>%
  as_tibble()

df_uk

Merge social distancing data

Identify London areas

nuts_london_inner <- c('UKI31','UKI32','UKI33','UKI34','UKI41',
                      'UKI42','UKI43','UKI44','UKI45')

nuts_london_outer <- c('UKI51','UKI52','UKI53','UKI54','UKI61',
                      'UKI62','UKI63','UKI71','UKI72','UKI73',
                      'UKI74','UKI75')

ut_london_inner <- c('E09000007','E09000001','E09000033','E09000013',
                    'E09000020','E09000032','E09000025','E09000012',
                    'E09000030','E09000014','E09000019','E09000023',
                    'E09000028','E09000022')

ut_london_outer <- c('E09000011','E09000004','E09000016','E09000002',
                    'E09000031','E09000026','E09000010','E09000006',
                    'E09000008','E09000029','E09000021','E09000024',
                    'E09000003','E09000005','E09000009','E09000017',
                    'E09000015','E09000018','E09000027')

df_uk = df_uk %>% 
  mutate(london = ifelse(ut_area %in% ut_london_inner, 'london_inner', 
                       ifelse(ut_area %in% ut_london_outer, 'london_outer',
                              'country'))) %>%
  mutate(london = as.factor(london))

df_uk_socdist = df_uk_socdist %>% 
  mutate(london = ifelse(nuts3 %in% nuts_london_inner, 'london_inner', 
                       ifelse(nuts3 %in% nuts_london_outer, 'london_outer',
                              'country'))) %>%
  mutate(london = as.factor(london))

Explore data

Plot prevalence over time

df_uk %>% ggplot(aes(x=time, y=rate_day)) + 
  geom_point(aes(col=ut_area, size=popdens)) + 
  geom_smooth(method="loess", se=T) + 
  theme(legend.position="none") +
  ggtitle("Overall prevalence over time")

pers <- c('pers_o', 'pers_c', 'pers_e', 'pers_a', 'pers_n')

for (i in pers){

gg <- df_uk %>% mutate(prev_tail = cut(.[[i]], 
                                       breaks = c(-Inf, quantile(.[[i]], 0.2), quantile(.[[i]], 0.8), Inf),
                                       labels = c('lower tail', 'center', 'upper tail'))) %>% 
  filter(prev_tail != 'center') %>%
  ggplot(aes(x=time, y=rate_day)) + 
  geom_point(aes(col=ut_area, size=popdens)) + 
  geom_smooth(method="loess", se=T) + 
  facet_wrap(~prev_tail) + 
  theme(legend.position="none") +
  ggtitle(i)

print(gg)
}

Plot social distancing over time

df_uk_socdist %>% ggplot(aes(x=time, y=socdist_single_tile)) + 
  geom_point(aes(col=nuts3, size=popdens)) + 
  geom_smooth(method="loess", se=T) + 
  theme(legend.position="none") +
  ggtitle("Overall social distancing over time")

pers <- c('pers_o', 'pers_c', 'pers_e', 'pers_a', 'pers_n')

for (i in pers){

gg <- df_uk_socdist %>% mutate(socdist_tail = cut(.[[i]], 
                                       breaks = c(-Inf, quantile(.[[i]], 0.2), quantile(.[[i]], 0.8), Inf),
                                       labels = c('lower tail', 'center', 'upper tail'))) %>% 
  filter(socdist_tail != 'center') %>%
  ggplot(aes(x=time, y=socdist_single_tile)) + 
  geom_point(aes(col=nuts3, size=popdens)) + 
  geom_smooth(method="loess", se=T) + 
  facet_wrap(~socdist_tail) + 
  theme(legend.position="none") +
  ggtitle(i)

print(gg)
}

Explore differences between london and the rest

df_uk %>% ggplot(aes(x=time, y=rate_day)) + 
  geom_point(aes(col=ut_area, size=popdens)) + 
  geom_smooth(method="loess", se=T) + 
  theme(legend.position="none") +
  facet_wrap(~london) +
  ggtitle("Overall prevalence over time")

NA
NA
NA

df_uk_socdist %>% ggplot(aes(x=time, y=socdist_single_tile)) + 
  geom_point(aes(col=nuts3, size=popdens)) + 
  geom_smooth(method="loess", se=T) + 
  theme(legend.position="none") +
  facet_wrap(~london) +
  ggtitle("Overall social distancing over time")

Control for weekend effect

df_uk_socdist %>% ggplot(aes(x=time, y=socdist_single_tile_clean)) + 
  geom_point(aes(col=nuts3, size=popdens)) + 
  geom_smooth(method="loess", se=T) + 
  theme(legend.position="none") +
  ggtitle("Overall social distancing over time")

pers <- c('pers_o', 'pers_c', 'pers_e', 'pers_a', 'pers_n')

for (i in pers){

gg <- df_uk_socdist %>% mutate(socdist_tail = cut(.[[i]], 
                                       breaks = c(-Inf, quantile(.[[i]], 0.2), quantile(.[[i]], 0.8), Inf),
                                       labels = c('lower tail', 'center', 'upper tail'))) %>% 
  filter(socdist_tail != 'center') %>%
  ggplot(aes(x=time, y=socdist_single_tile_clean)) + 
  geom_point(aes(col=nuts3, size=popdens)) + 
  geom_smooth(method="loess", se=T) + 
  facet_wrap(~socdist_tail) + 
  theme(legend.position="none") +
  ggtitle(i)

print(gg)
}

df_uk_socdist <- df_uk_socdist %>% mutate(socdist_single_tile = socdist_single_tile_clean) %>% 
  select(-loess, -socdist_single_tile_clean)

Correlations

df_uk %>% group_by(ut_area) %>% 
  summarize_if(is.numeric, mean, na.rm=T) %>% 
  select(-ut_area, -time) %>% 
  cor(use = 'pairwise.complete') %>% round(3)

              rate_day pers_o pers_a pers_n pers_c pers_e airport_dist  males popdens manufacturing tourism health
rate_day         1.000  0.249 -0.256  0.080 -0.425  0.199       -0.452  0.035   0.420        -0.197  -0.280  0.133
pers_o           0.249  1.000 -0.604 -0.143 -0.636  0.734       -0.247  0.437   0.806        -0.532   0.116  0.089
pers_a          -0.256 -0.604  1.000 -0.197  0.643 -0.412        0.285 -0.460  -0.745         0.496   0.069 -0.001
pers_n           0.080 -0.143 -0.197  1.000 -0.353 -0.497        0.059  0.096   0.043         0.359  -0.244  0.138
pers_c          -0.425 -0.636  0.643 -0.353  1.000 -0.293        0.327 -0.524  -0.727         0.363   0.227 -0.286
pers_e           0.199  0.734 -0.412 -0.497 -0.293  1.000       -0.262  0.197   0.617        -0.595   0.098  0.026
airport_dist    -0.452 -0.247  0.285  0.059  0.327 -0.262        1.000 -0.124  -0.381         0.300   0.448  0.150
males            0.035  0.437 -0.460  0.096 -0.524  0.197       -0.124  1.000   0.510        -0.259  -0.094 -0.012
popdens          0.420  0.806 -0.745  0.043 -0.727  0.617       -0.381  0.510   1.000        -0.545  -0.099  0.112
manufacturing   -0.197 -0.532  0.496  0.359  0.363 -0.595        0.300 -0.259  -0.545         1.000  -0.088 -0.023
tourism         -0.280  0.116  0.069 -0.244  0.227  0.098        0.448 -0.094  -0.099        -0.088   1.000  0.049
health           0.133  0.089 -0.001  0.138 -0.286  0.026        0.150 -0.012   0.112        -0.023   0.049  1.000
academic         0.219  0.707 -0.502 -0.447 -0.298  0.730       -0.369  0.267   0.602        -0.714   0.056 -0.228
medinc           0.235  0.533 -0.533 -0.335 -0.300  0.569       -0.379  0.367   0.631        -0.532  -0.062 -0.243
medage          -0.345 -0.508  0.582 -0.167  0.780 -0.337        0.447 -0.620  -0.705         0.472   0.435 -0.153
conservative    -0.260 -0.843  0.511  0.366  0.485 -0.767        0.396 -0.293  -0.686         0.657  -0.023  0.043
              academic medinc medage conservative
rate_day         0.219  0.235 -0.345       -0.260
pers_o           0.707  0.533 -0.508       -0.843
pers_a          -0.502 -0.533  0.582        0.511
pers_n          -0.447 -0.335 -0.167        0.366
pers_c          -0.298 -0.300  0.780        0.485
pers_e           0.730  0.569 -0.337       -0.767
airport_dist    -0.369 -0.379  0.447        0.396
males            0.267  0.367 -0.620       -0.293
popdens          0.602  0.631 -0.705       -0.686
manufacturing   -0.714 -0.532  0.472        0.657
tourism          0.056 -0.062  0.435       -0.023
health          -0.228 -0.243 -0.153        0.043
academic         1.000  0.733 -0.382       -0.888
medinc           0.733  1.000 -0.481       -0.661
medage          -0.382 -0.481  1.000        0.486
conservative    -0.888 -0.661  0.486        1.000

df_uk_socdist %>% group_by(nuts3) %>% 
  summarize_if(is.numeric, mean, na.rm=T) %>% 
  select(-nuts3, -time) %>% 
  cor(use = 'pairwise.complete') %>% round(3)

                    socdist_tiles socdist_single_tile pers_o pers_e pers_a pers_n pers_c airport_dist  males popdens
socdist_tiles               1.000              -0.467 -0.632 -0.735  0.546  0.337  0.397        0.530 -0.296  -0.711
socdist_single_tile        -0.467               1.000  0.189  0.235 -0.276  0.091 -0.211       -0.270 -0.004   0.376
pers_o                     -0.632               0.189  1.000  0.717 -0.604 -0.128 -0.655       -0.231  0.522   0.805
pers_e                     -0.735               0.235  0.717  1.000 -0.460 -0.451 -0.374       -0.319  0.246   0.644
pers_a                      0.546              -0.276 -0.604 -0.460  1.000 -0.203  0.659        0.328 -0.568  -0.748
pers_n                      0.337               0.091 -0.128 -0.451 -0.203  1.000 -0.324        0.072  0.190   0.020
pers_c                      0.397              -0.211 -0.655 -0.374  0.659 -0.324  1.000        0.326 -0.630  -0.737
airport_dist                0.530              -0.270 -0.231 -0.319  0.328  0.072  0.326        1.000 -0.155  -0.377
males                      -0.296              -0.004  0.522  0.246 -0.568  0.190 -0.630       -0.155  1.000   0.613
popdens                    -0.711               0.376  0.805  0.644 -0.748  0.020 -0.737       -0.377  0.613   1.000
manufacturing               0.676              -0.349 -0.497 -0.607  0.430  0.351  0.310        0.272 -0.179  -0.492
tourism                     0.225              -0.055  0.083 -0.011  0.129 -0.137  0.194        0.492 -0.168  -0.130
health                      0.132              -0.034  0.087  0.045  0.004  0.155 -0.275        0.180  0.023   0.166
academic                   -0.802               0.198  0.713  0.742 -0.491 -0.432 -0.324       -0.379  0.263   0.574
medinc                     -0.767               0.320  0.552  0.582 -0.620 -0.247 -0.367       -0.389  0.440   0.652
medage                      0.569              -0.218 -0.512 -0.402  0.631 -0.181  0.790        0.469 -0.672  -0.714
conservative                0.751              -0.172 -0.837 -0.761  0.514  0.338  0.522        0.389 -0.350  -0.678
rate_day                   -0.705               0.361  0.591  0.563 -0.549 -0.125 -0.457       -0.397  0.267   0.695
                    manufacturing tourism health academic medinc medage conservative rate_day
socdist_tiles               0.676   0.225  0.132   -0.802 -0.767  0.569        0.751   -0.705
socdist_single_tile        -0.349  -0.055 -0.034    0.198  0.320 -0.218       -0.172    0.361
pers_o                     -0.497   0.083  0.087    0.713  0.552 -0.512       -0.837    0.591
pers_e                     -0.607  -0.011  0.045    0.742  0.582 -0.402       -0.761    0.563
pers_a                      0.430   0.129  0.004   -0.491 -0.620  0.631        0.514   -0.549
pers_n                      0.351  -0.137  0.155   -0.432 -0.247 -0.181        0.338   -0.125
pers_c                      0.310   0.194 -0.275   -0.324 -0.367  0.790        0.522   -0.457
airport_dist                0.272   0.492  0.180   -0.379 -0.389  0.469        0.389   -0.397
males                      -0.179  -0.168  0.023    0.263  0.440 -0.672       -0.350    0.267
popdens                    -0.492  -0.130  0.166    0.574  0.652 -0.714       -0.678    0.695
manufacturing               1.000  -0.052 -0.067   -0.640 -0.462  0.413        0.627   -0.403
tourism                    -0.052   1.000  0.024    0.005 -0.138  0.485        0.012   -0.087
health                     -0.067   0.024  1.000   -0.202 -0.261 -0.154        0.010   -0.012
academic                   -0.640   0.005 -0.202    1.000  0.724 -0.368       -0.891    0.624
medinc                     -0.462  -0.138 -0.261    0.724  1.000 -0.479       -0.653    0.595
medage                      0.413   0.485 -0.154   -0.368 -0.479  1.000        0.499   -0.475
conservative                0.627   0.012  0.010   -0.891 -0.653  0.499        1.000   -0.610
rate_day                   -0.403  -0.087 -0.012    0.624  0.595 -0.475       -0.610    1.000

Modelling

Prepare functions

# function calculates all relevant models
run_models <- function(y, lvl1_x, lvl2_x, lvl2_id, data, ctrls=F){

  # subset data
  data = data %>% 
    dplyr::select(all_of(y), all_of(lvl1_x), all_of(lvl2_x), all_of(lvl2_id), 
                  popdens, rate_day, all_of(y))
  data = data %>% 
    dplyr::rename(y = all_of(y),
           lvl1_x = all_of(lvl1_x),
           lvl2_x = all_of(lvl2_x),
           lvl2_id = all_of(lvl2_id)
           )
  
  # configure optimization procedure
  ctrl_config <- lmeControl(opt = 'optim', maxIter = 100, msMaxIter = 100)

  # baseline
  baseline <- lme(fixed = y ~ 1, random = ~ 1 | lvl2_id, 
                    data = data,
                    correlation = corAR1(),
                  control = ctrl_config,
                  method = 'ML')

  # random intercept fixed slope
  random_intercept <- lme(fixed = y ~ lvl1_x + lvl2_x, 
                          random = ~ 1 | lvl2_id,
                            data = data,
                            correlation = corAR1(),
                  control = ctrl_config,
                  method = 'ML')

  # random intercept random slope
  random_slope <- lme(fixed = y ~ lvl1_x + lvl2_x, 
                      random = ~ lvl1_x | lvl2_id, 
                        data = data,
                        correlation = corAR1(),
                  control = ctrl_config,
                  method = 'ML')

  # cross level interaction
  interaction <- lme(fixed = y ~ lvl1_x * lvl2_x, 
                     random = ~ lvl1_x | lvl2_id, 
                       data = data,
                       correlation = corAR1(),
                  control = ctrl_config,
                  method = 'ML')
  
  # create list with results
  results <- list('baseline' = baseline, 
                  "random_intercept" = random_intercept, 
                  "random_slope" = random_slope,
                  "interaction" = interaction)
  
  
  if (ctrls == 'dem' | ctrls == 'prev'){
    
    # random intercept random slope
    random_slope_ctrl_dem <- lme(fixed = y ~ lvl1_x + lvl2_x + popdens,
                              random = ~ lvl1_x | lvl2_id, 
                          data = data,
                          correlation = corAR1(),
                    control = ctrl_config,
                  method = 'ML')
  
    # cross level interaction
    interaction_ctrl_main_dem <- lme(fixed = y ~ lvl1_x * lvl2_x + popdens,
                             random = ~ lvl1_x | lvl2_id, 
                         data = data,
                         correlation = corAR1(),
                    control = ctrl_config,
                  method = 'ML')
  
    # cross level interaction
    interaction_ctrl_int_dem <- lme(fixed = y ~ lvl1_x * lvl2_x + lvl1_x * popdens,
                             random = ~ lvl1_x | lvl2_id, 
                         data = data,
                         correlation = corAR1(),
                    control = ctrl_config,
                  method = 'ML')        
    
    # create list with results
    results <- list('baseline' = baseline, 
                    "random_intercept" = random_intercept, 
                    "random_slope" = random_slope,
                    "interaction" = interaction,
                    "random_slope_ctrl_dem" = random_slope_ctrl_dem,
                    "interaction_ctrl_main_dem" = interaction_ctrl_main_dem,
                    "interaction_ctrl_int_dem" = interaction_ctrl_int_dem)
  }
  
  if (ctrls == 'prev'){
  
    # random intercept random slope
    random_slope_ctrl_prev <- lme(fixed = y ~ lvl1_x + lvl2_x + popdens + rate_day,
                              random = ~ lvl1_x + rate_day | lvl2_id, 
                          data = data,
                          correlation = corAR1(),
                          control = ctrl_config,
                  method = 'ML')  
    
        # cross level interaction
    interaction_ctrl_main_prev <- lme(fixed = y ~ lvl1_x * lvl2_x + popdens + rate_day,
                             random = ~ lvl1_x | lvl2_id, 
                         data = data,
                         correlation = corAR1(),
                    control = ctrl_config,
                  method = 'ML')
  
  
    # cross level interaction
    interaction_ctrl_int_prev<- lme(fixed = y ~ lvl1_x * lvl2_x + lvl1_x * popdens + rate_day,
                             random = ~ lvl1_x + rate_day | lvl2_id, 
                         data = data,
                         correlation = corAR1(),
                          control = ctrl_config,
                  method = 'ML')
  
    # create list with results
    results <- list('baseline' = baseline, 
                    "random_intercept" = random_intercept, 
                    "random_slope" = random_slope,
                    "interaction" = interaction,
                    "random_slope_ctrl_dem" = random_slope_ctrl_dem,
                    "interaction_ctrl_main_dem" = interaction_ctrl_main_dem,
                    "interaction_ctrl_int_dem" = interaction_ctrl_int_dem,                    
                    "random_slope_ctrl_prev" = random_slope_ctrl_prev,
                    "interaction_ctrl_main_prev" = interaction_ctrl_main_prev,
                    "interaction_ctrl_int_prev" = interaction_ctrl_int_prev)
  }
  
  if(ctrls == 'exp'){
    # random intercept random slope
  random_slope_exp <- lme(fixed = y ~ (lvl1_x + I(lvl1_x^2)) + lvl2_x, 
                      random = ~ (lvl1_x + I(lvl1_x^2)) | lvl2_id, 
                        data = data,
                        correlation = corAR1(),
                  control = ctrl_config,
                  method = 'ML')

  # cross level interaction
  interaction_exp <- lme(fixed = y ~ (lvl1_x + I(lvl1_x^2)) * lvl2_x, 
                     random = ~ (lvl1_x + I(lvl1_x^2)) | lvl2_id, 
                       data = data,
                       correlation = corAR1(),
                  control = ctrl_config,
                  method = 'ML')  
  
  
  # create list with results
  results <- list('baseline' = baseline, 
                  "random_intercept" = random_intercept, 
                  "random_slope" = random_slope,
                  "interaction" = interaction,                  
                  "random_slope_exp" = random_slope_exp,
                  "interaction_exp" = interaction_exp)
  }
  
  return(results)
        
}

# extracts table with coefficients and tests statistics
extract_results <- function(models) {
  
  models_summary <- models %>% 
  map(summary) %>% 
  map("tTable") %>% 
  map(as.data.frame) %>% 
  map(round, 10) 
  # %>% map(~ .[str_detect(rownames(.), 'Inter|lvl'),])
  
  return(models_summary)
  
}


# calculates comparison of all models in model list
compare_models <- function(models) {

  mdl_names <- models %>% names()
  
  str = ''
  for (i in mdl_names){
    
    mdl_str <- paste('models$', i, sep = '')
    
    if(i == 'baseline'){
      str <- mdl_str
    }else{
    str <- paste(str, mdl_str, sep=', ')
    }
  }
  
  anova_str <- paste0('anova(', str, ')')
  mdl_comp <- eval(parse(text=anova_str))
  rownames(mdl_comp) = mdl_names
  return(mdl_comp)
}

Remove London Data

# df_uk <- df_uk %>% filter(london == 'country')
# df_uk_socdist <- df_uk_socdist %>% filter(london == 'country')

Rescale Data

Adjust timeframes

Predict prevalence

prevalence ~ openness

models_o_covid <-run_models(y = 'rate_day', 
                         lvl1_x = 'time', 
                         lvl2_x = 'pers_o', 
                         lvl2_id = 'ut_area', 
                         data = df_uk_scaled_prev,
                         ctrls = 'dem')

extract_results(models_o_covid)

compare_models(models_o_covid)

prevalence ~ conscientiousness

models_c_covid <-run_models(y = 'rate_day', 
                         lvl1_x = 'time', 
                         lvl2_x = 'pers_c', 
                         lvl2_id = 'ut_area', 
                         data = df_uk_scaled_prev,
                         ctrls = 'dem')

extract_results(models_c_covid)

compare_models(models_c_covid)

prevalence ~ extraversion

models_e_covid <-run_models(y = 'rate_day', 
                         lvl1_x = 'time', 
                         lvl2_x = 'pers_e', 
                         lvl2_id = 'ut_area', 
                         data = df_uk_scaled_prev,
                         ctrls = 'dem')

extract_results(models_e_covid)

compare_models(models_e_covid)

prevalence ~ agreeableness

models_a_covid <-run_models(y = 'rate_day', 
                         lvl1_x = 'time', 
                         lvl2_x = 'pers_a', 
                         lvl2_id = 'ut_area', 
                         data = df_uk_scaled_prev,
                         ctrls = 'dem')

extract_results(models_a_covid)

compare_models(models_a_covid)

prevalence ~ neuroticism

models_n_covid <-run_models(y = 'rate_day', 
                         lvl1_x = 'time', 
                         lvl2_x = 'pers_n', 
                         lvl2_id = 'ut_area', 
                         data = df_uk_scaled_prev,
                         ctrls = 'dem')

extract_results(models_n_covid)

compare_models(models_n_covid)

Predict social distancing

social distancing ~ openness

models_o_socdist <-run_models(y = 'socdist_single_tile', 
                         lvl1_x = 'time', 
                         lvl2_x = 'pers_o', 
                         lvl2_id = 'nuts3', 
                         data = df_uk_socdist_scaled,
                         ctrls = 'prev')

extract_results(models_o_socdist)

compare_models(models_o_socdist)

social distancing ~ conscientiousness

models_c_socdist <-run_models(y = 'socdist_single_tile', 
                         lvl1_x = 'time', 
                         lvl2_x = 'pers_c', 
                         lvl2_id = 'nuts3', 
                         data = df_uk_socdist_scaled,
                         ctrls = 'prev')

extract_results(models_c_socdist)

compare_models(models_c_socdist)

social distancing ~ extraversion

models_e_socdist <-run_models(y = 'socdist_single_tile', 
                         lvl1_x = 'time', 
                         lvl2_x = 'pers_e', 
                         lvl2_id = 'nuts3', 
                         data = df_uk_socdist_scaled,
                         ctrls = 'prev')

extract_results(models_e_socdist)

compare_models(models_e_socdist)

social distancing ~ agreeableness

models_a_socdist <-run_models(y = 'socdist_single_tile', 
                         lvl1_x = 'time', 
                         lvl2_x = 'pers_a', 
                         lvl2_id = 'nuts3', 
                         data = df_uk_socdist_scaled,
                         ctrls = 'prev')

extract_results(models_a_socdist)

compare_models(models_a_socdist)

social distancing ~ neuroticism

models_n_socdist <-run_models(y = 'socdist_single_tile', 
                         lvl1_x = 'time', 
                         lvl2_x = 'pers_n', 
                         lvl2_id = 'nuts3', 
                         data = df_uk_socdist_scaled,
                         ctrls = 'prev')

extract_results(models_n_socdist)

compare_models(models_n_socdist)

Explore quadratic trends

prevalence ~ openness

models_o_covid_exp <-run_models(y = 'rate_day',
                         lvl1_x = 'time',
                         lvl2_x = 'pers_o',
                         lvl2_id = 'ut_area',
                         data = df_uk_scaled_prev,
                         ctrls = 'exp')

extract_results(models_o_covid_exp)

compare_models(models_o_covid_exp)

prevalence ~ conscientiousness

models_c_covid_exp <-run_models(y = 'rate_day', 
                         lvl1_x = 'time', 
                         lvl2_x = 'pers_c', 
                         lvl2_id = 'ut_area', 
                         data = df_uk_scaled_prev,
                         ctrls = 'exp')

extract_results(models_c_covid_exp)

compare_models(models_c_covid_exp)

prevalence ~ extraversion

models_e_covid_exp <-run_models(y = 'rate_day', 
                         lvl1_x = 'time', 
                         lvl2_x = 'pers_e', 
                         lvl2_id = 'ut_area', 
                         data = df_uk_scaled_prev,
                         ctrls = 'exp')

extract_results(models_e_covid_exp)

compare_models(models_e_covid_exp)

prevalence ~ agreeableness

models_a_covid_exp <-run_models(y = 'rate_day', 
                         lvl1_x = 'time', 
                         lvl2_x = 'pers_a', 
                         lvl2_id = 'ut_area', 
                         data = df_uk_scaled_prev,
                         ctrls = 'exp')

extract_results(models_a_covid_exp)

compare_models(models_a_covid_exp)

prevalence ~ neuroticism

models_n_covid_exp <-run_models(y = 'rate_day', 
                         lvl1_x = 'time', 
                         lvl2_x = 'pers_n', 
                         lvl2_id = 'ut_area', 
                         data = df_uk_scaled_prev,
                         ctrls = 'exp')

extract_results(models_n_covid_exp)

compare_models(models_n_covid_exp)

Create overview table

Define function to create overview tables

summary_table <- function(models, dv_name, prev=F){

  temp_df_ctrl_main <- NULL
  temp_df_ctrl_int <- NULL
  temp_df_ctrl_int_prev <- NULL
  
  for (i in models){
    results <- i %>% extract_results()
    
    results_ctrl_main <- results$interaction_ctrl_main_dem['lvl1_x:lvl2_x',]
    temp_df_ctrl_main <- temp_df_ctrl_main %>% rbind(results_ctrl_main)
    
    results_ctrl_int <- results$interaction_ctrl_int_dem['lvl1_x:lvl2_x',]
    temp_df_ctrl_int <- temp_df_ctrl_int %>% rbind(results_ctrl_int)
    
    if(prev){
      results_ctrl_int_prev <- results$interaction_ctrl_int_prev['lvl1_x:lvl2_x',]
      temp_df_ctrl_int_prev <- temp_df_ctrl_int_prev %>% rbind(results_ctrl_int_prev)
    }
        
  }
  
  names_ctrl_main <- paste0(dv_name, '~', c('o', 'c', 'e', 'a', 'n'), '*time', '_crtl_popdens')
  rownames(temp_df_ctrl_main) <- names_ctrl_main

  names_ctrl_int <- paste0(dv_name, '~', c('o', 'c', 'e', 'a', 'n'), '*time', '_crtl_popdens*time')
  rownames(temp_df_ctrl_int) <- names_ctrl_int

  if(prev){
    names_ctrl_int_prev <- paste0(dv_name, '~', c('o', 'c', 'e', 'a', 'n'), '*time', '_crtl_popdens*time_prev')
    rownames(temp_df_ctrl_int_prev) <- names_ctrl_int_prev
    
    sum_tab <- rbind(temp_df_ctrl_main, temp_df_ctrl_int, temp_df_ctrl_int_prev) %>% round(4)
  }else{
    sum_tab <- rbind(temp_df_ctrl_main, temp_df_ctrl_int) %>% round(4)
  }


  
  return(sum_tab)

} 

Create overview tables

# prevalence
models_prev <- list(models_o_covid, 
                    models_c_covid, 
                    models_e_covid, 
                    models_a_covid, 
                    models_n_covid)

sum_tab_prev <- summary_table(models_prev, dv_name = 'prev')

write.table(sum_tab_prev, quote=F)

# social distancing
models_socdist <- list(models_o_socdist, 
                       models_c_socdist, 
                       models_e_socdist, 
                       models_a_socdist, 
                       models_n_socdist)

sum_tab_socdist <- summary_table(models_socdist, dv_name = 'socdist', prev=T)

write.table(sum_tab_socdist, quote=F)

Conditional random forest analysis

Extract slopes prevalence

# slope prevalence
df_uk_slope_prev <- df_uk_scaled_prev %>% split(.$ut_area) %>% 
  map(~ lm(rate_day ~ time, data = .)) %>%
  map(coef) %>% 
  map_dbl('time') %>% 
  as.data.frame() %>% 
  rownames_to_column('ut_area') %>% 
  rename(slope_prev = '.')

# merge with control variables 
df_uk_slope_prev <- df_uk_scaled_prev %>% select(-time, -rate_day) %>%
  distinct() %>% 
  inner_join(df_uk_slope_prev, by = 'ut_area') %>%
  drop_na()

head(df_uk_slope_prev)

Extract slopes social distancing

# slope socdist
df_uk_slope_socdist <- df_uk_socdist_scaled %>% split(.$nuts3) %>%
  map(~ lm(socdist_single_tile ~ time, data = .)) %>%
  map(coef) %>%
  map_dbl('time') %>%
  as.data.frame() %>%
  rownames_to_column('nuts3') %>%
  rename(slope_socdist = '.')

# merge with control variables 
df_uk_slope_socdist <- df_uk_socdist_scaled %>% 
  select(-time, -date, -socdist_tiles, -socdist_single_tile) %>%
  distinct() %>%
  inner_join(df_uk_slope_socdist, by = 'nuts3') %>%
  drop_na()

head(df_uk_slope_socdist)

Explore distribution of slopes

df_uk_slope_prev %>% ggplot(aes(slope_prev)) + geom_histogram(bins = 100)

df_uk_slope_socdist %>% ggplot(aes(slope_socdist)) + geom_histogram(bins = 100)

df_uk_slope_prev

CRF prevalence ~ openness

ctrls <- cforest_unbiased(ntree=500, mtry=5)

crf_o_fit_prev <- cforest(slope_prev ~ pers_o + airport_dist + males +
                          popdens + manufacturing + tourism +
                          health + academic + medinc + medage + conservative, 
                         df_uk_slope_prev[-1], 
                         controls = ctrls)

crf_o_varimp_prev <- varimp(crf_o_fit_prev, nperm = 1)
crf_o_varimp_cond_prev <- varimp(crf_o_fit_prev, conditional = T, nperm = 1)

crf_o_varimp_prev %>% as.data.frame() %>% rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))

crf_o_varimp_cond_prev %>% as.data.frame() %>% rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') +
  theme(axis.text.x = element_text(angle = 90))

CRF prevalence ~ conscientiousness

ctrls <- cforest_unbiased(ntree=500, mtry=5)

crf_c_fit_prev <- cforest(slope_prev ~ pers_c + airport_dist + males +
                          popdens + manufacturing + tourism +
                          health + academic + medinc + medage + conservative, 
                         df_uk_slope_prev[-1], 
                         controls = ctrls)

crf_c_varimp_prev <- varimp(crf_c_fit_prev, nperm = 1)
crf_c_varimp_cond_prev <- varimp(crf_c_fit_prev, conditional = T, nperm = 1)

crf_c_varimp_prev %>% as.data.frame() %>% rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))

crf_c_varimp_cond_prev %>% as.data.frame() %>% rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))

CRF prevalence ~ extraversion

ctrls <- cforest_unbiased(ntree=500, mtry=5)

crf_e_fit_prev <- cforest(slope_prev ~ pers_e + airport_dist + males +
                          popdens + manufacturing + tourism +
                          health + academic + medinc + medage + conservative, 
                         df_uk_slope_prev[-1], 
                         controls = ctrls)

crf_e_varimp_prev <- varimp(crf_e_fit_prev, nperm = 1)
crf_e_varimp_cond_prev <- varimp(crf_e_fit_prev, conditional = T, nperm = 1)

crf_e_varimp_prev %>% as.data.frame() %>% rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))

crf_e_varimp_cond_prev %>% as.data.frame() %>% rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))

CRF prevalence ~ agreeableness

ctrls <- cforest_unbiased(ntree=500, mtry=5)

crf_a_fit_prev <- cforest(slope_prev ~ pers_a + airport_dist + males +
                          popdens + manufacturing + tourism +
                          health + academic + medinc + medage + conservative, 
                         df_uk_slope_prev[-1], 
                         controls = ctrls)

crf_a_varimp_prev <- varimp(crf_a_fit_prev, nperm = 1)
crf_a_varimp_cond_prev <- varimp(crf_a_fit_prev, conditional = T, nperm = 1)

crf_a_varimp_prev %>% as.data.frame() %>% rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))

crf_a_varimp_cond_prev %>% as.data.frame() %>% rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))

CRF prevalence ~ neuroticism

ctrls <- cforest_unbiased(ntree=500, mtry=5)

crf_n_fit_prev <- cforest(slope_prev ~ pers_n + airport_dist + males +
                          popdens + manufacturing + tourism +
                          health + academic + medinc + medage + conservative, 
                         df_uk_slope_prev[-1], 
                         controls = ctrls)

crf_n_varimp_prev <- varimp(crf_n_fit_prev, nperm = 1)
crf_n_varimp_cond_prev <- varimp(crf_n_fit_prev, conditional = T, nperm = 1)

crf_n_varimp_prev %>% as.data.frame() %>% rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))

crf_n_varimp_cond_prev %>% as.data.frame() %>% rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))

CRF social distancing ~ openness

ctrls <- cforest_unbiased(ntree=500, mtry=5)

crf_o_fit_socdist <- cforest(slope_socdist ~ pers_o + airport_dist + males +
                          popdens + manufacturing + tourism +
                          health + academic + medinc + medage + conservative, 
                         df_uk_slope_socdist[-1], 
                         controls = ctrls)

crf_o_varimp_socdist <- varimp(crf_o_fit_socdist, nperm = 1)
crf_o_varimp_cond_socdist <- varimp(crf_o_fit_socdist, conditional = T, nperm = 1)

crf_o_varimp_socdist %>% as.data.frame() %>% rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))

crf_o_varimp_cond_socdist %>% as.data.frame() %>% rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))

CRF social distancing ~ conscientiousness

ctrls <- cforest_unbiased(ntree=500, mtry=5)

crf_c_fit_socdist <- cforest(slope_socdist ~ pers_c + airport_dist + males +
                          popdens + manufacturing + tourism +
                          health + academic + medinc + medage + conservative, 
                         df_uk_slope_socdist[-1], 
                         controls = ctrls)

crf_c_varimp_socdist <- varimp(crf_c_fit_socdist, nperm = 1)
crf_c_varimp_cond_socdist <- varimp(crf_c_fit_socdist, conditional = T, nperm = 1)

crf_c_varimp_socdist %>% as.data.frame() %>% rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))

crf_c_varimp_cond_socdist %>% as.data.frame() %>% rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))

CRF social distancing ~ extraversion

ctrls <- cforest_unbiased(ntree=500, mtry=5)

crf_e_fit_socdist <- cforest(slope_socdist ~ pers_e + airport_dist + males +
                          popdens + manufacturing + tourism +
                          health + academic + medinc + medage + conservative, 
                         df_uk_slope_socdist[-1], 
                         controls = ctrls)

crf_e_varimp_socdist <- varimp(crf_e_fit_socdist, nperm = 1)
crf_e_varimp_cond_socdist <- varimp(crf_e_fit_socdist, conditional = T, nperm = 1)

crf_e_varimp_socdist %>% as.data.frame() %>% rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))

crf_e_varimp_cond_socdist %>% as.data.frame() %>% rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))

CRF social distancing ~ agreeableness

ctrls <- cforest_unbiased(ntree=500, mtry=5)

crf_a_fit_socdist <- cforest(slope_socdist ~ pers_a + airport_dist + males +
                          popdens + manufacturing + tourism +
                          health + academic + medinc + medage + conservative, 
                         df_uk_slope_socdist[-1], 
                         controls = ctrls)

crf_a_varimp_socdist <- varimp(crf_a_fit_socdist, nperm = 1)
crf_a_varimp_cond_socdist <- varimp(crf_a_fit_socdist, conditional = T, nperm = 1)

crf_a_varimp_socdist %>% as.data.frame() %>% rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))

crf_a_varimp_cond_socdist %>% as.data.frame() %>% rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))

CRF social distancing ~ neuroticism

ctrls <- cforest_unbiased(ntree=500, mtry=5)

crf_n_fit_socdist <- cforest(slope_socdist ~ pers_n + airport_dist + males +
                          popdens + manufacturing + tourism +
                          health + academic + medinc + medage + conservative, 
                         df_uk_slope_socdist[-1], 
                         controls = ctrls)

crf_n_varimp_socdist <- varimp(crf_n_fit_socdist, nperm = 1)
crf_n_varimp_cond_socdist <- varimp(crf_n_fit_socdist, conditional = T, nperm = 1)

crf_n_varimp_socdist %>% as.data.frame() %>% rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))

crf_n_varimp_cond_socdist %>% as.data.frame() %>% rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))

Change point analysis

Preparation

# keep only counties with full data
ut_area_complete <- df_uk_scaled %>% 
  group_by(ut_area) %>%
  summarize(n = n()) %>%
  filter(n==max(.$n)) %>% 
  .$ut_area

Prevalence

df_uk_cpt_prev %>% select(cpt_day_prev) %>% map(hist)

$cpt_day_prev
$breaks
[1] 20 25 30 35 40 45 50

$counts
[1] 107   0   5  14  21   2

$density
[1] 0.143624161 0.000000000 0.006711409 0.018791946 0.028187919 0.002684564

$mids
[1] 22.5 27.5 32.5 37.5 42.5 47.5

$xname
[1] ".x[[i]]"

$equidist
[1] TRUE

attr(,"class")
[1] "histogram"

Social distancing

for(i in head(df_uk_socdist_cpt_results,5)){
  plot(i)
}

Predicting change points

Linear models predicting change points (no controls)

lm_cpr_prev_pers <- lm(cpt_day_prev ~ pers_o + pers_c + pers_e + pers_a + pers_n, 
                         data = df_uk_cpt_prev_socdist)
lm_cpr_prev_pers %>% summary()


lm_cpt_socdist_pers <- lm(cpt_day_socdist ~ pers_o + pers_c + pers_e + pers_a + pers_n, 
                            data = df_uk_cpt_prev_socdist)
lm_cpt_socdist_pers %>% summary()

Linear models predicting change points with controls

df_uk_cpt_prev_socdist

lm_cpt_prev_pers <- lm(cpt_day_prev ~ pers_o + pers_c + pers_e + pers_a + pers_n + 
                         women + academics + hospital_beds + gdp + manufact +
                          airport + age + popdens,
                         data = df_uk_cpt_prev_socdist)
lm_cpt_prev_pers %>% summary()

lm_cpt_socdist_pers <- lm(cpt_day_socdist ~ pers_o + pers_c + pers_e + pers_a + pers_n + 
                            women + academics + hospital_beds + gdp + manufact +
                            airport + age + popdens,
                            data = df_uk_cpt_prev_socdist)
lm_cpt_socdist_pers %>% summary()